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Scalar-tensor theories are a compelling alternative to general relativity and one of the most ac- 
cepted extensions of Einstein's theory. Black holes in these theories have no hair, but could grow 
"wigs" supported by time-dependent boundary conditions or spatial gradients. Time-dependent or 
spatially varying fields lead in general to nontrivial black hole dynamics, with potentially interesting 
experimental consequences. We carry out a numerical investigation of the dynamics of single and 
binary black holes in the presence of scalar fields. In particular we study gravitational and scalar 
radiation from black- hole binaries in a constant scalar-field gradient, and we compare our numerical 
findings to analytical models. 
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I. INTRODUCTION 

Scalar fields are ubiquitous in physics, either as a proxy 
for more complex interactions or as fundamental quan- 
tities in their own right. For example, one of the best 
studied modifications to general relativity is scalar-tensor 
gravity, in which space-time curvature couples to scalar 
fields (which are sufficiently light to be relevant for astro- 
physics and/or cosmology). The historical development 
of this theory goes back to the 1940s, and involves several 
research groups with different views on the physical inter- 
pretation of the scalar degree of freedom [IHS] . In recent 
times, interest in scalar-tensor gravity has been driven 
by theoretical attempts to unify gravity with quantum 
mechanics at high energies and solve the cosmological 
constant and hierarchy problems, as well as observations 
of the cosmic microwave background and the highly- 
anticipated future observations of gravitational waves [4]- 

The simplest version of scalar-tensor gravity is Brans- 
Dicke theory JHP], which consists of a single massless 
scalar (f), whose coupling to curvature is controlled by a 
dimensionless parameter wbd- Generalizations include 
varying scalar-curvature couplings uj{4>) and a scalar po- 
tential ¥{(/)) ( "Bergmann- Wagoner" theories [IHllIl]) as 
well as the possibility of multiple interacting scalar fields: 
see e.g. [H [T^l [13] for comprehensive treatments of the 
subject. These generalizations (commonly referred to as 
"scalar-tensor theories" ) are compelling due to their sim- 
plicity, but (perhaps as a consequence) they are also very 



|berti@phy.olemiss.edu| 



well constrained observationally. The classic book by 
Will [141 presents an overview of the subject, and com- 
prehensive reviews of the state of the art in experimental 
tests of gravitational theories can be found in ^E\ [1^1 . 

In this paper we shall focus, for simplicity, on scalar- 
tensor theories involving a single scalar field. In general 
relativity, because of the conservation of total momentum 
for isolated systems, gravitational radiation is quadrupo- 
lar in nature. A possible smoking gun of scalar-tensor 
gravity is the existence of dipole radiation, essentially 
due to violations of the equivalence principle. Solar Sys- 
tem experiments and observations of binary pulsar sys- 
tems place strong constraints on the coupling functions of 
the theory [HI [TH [TSl [T7H27] . and other interesting con- 
straints may come from the direct observation of grav- 
itational radiation from binary systems in the near fu- 
ture [m [25H55] . Despite all of these observational con- 
straints, striking and potentially observable astrophysical 
phenomena are still possible in these theories. Such phe- 
nomena include superradiant instabilities (see e.g. [5S1 - 
I39j for discussions of this phenomenon in the context of 
the "string axiverse" scenario [40:) and the related possi- 
bility of floating orbits around rotating black holes (BHs) 



A. Classical no-hair theorems 

Theoretical studies impose remarkable constraints and 
limitations on scalar-tensor theories. First of all, the fa- 
mous BH no-scalar-hair theorems first proved by Hawk- 
ing g2], Thorne and Dykla [43], and Chase [ll] state 
that stationary BH solutions in Brans-Dicke theory are 



the same as those in general relativity. These results 
have been generalized and expanded upon by many au- 
thors. For example, an extension to multiple scalars has 
been established by Heusler "iS], and an extension to 
Bergmann- Wagoner and f{R) theories has been estab- 
lished by Sotiriou and Faraoni [46 . The no-scalar- hair 
theorems have also been confirmed by numerical stud- 
ies of gravitational collapse ;47-n52| . More generally, it 
has been observed that the Kerr metric is a solution in 
a wide class of gravity theories [53]. However, the fact 
that stationary vacuum solutions of scalar-tensor theo- 
ries agree with those of general relativity does not mean 
that the dynamics of BHs in these theories must be the 
same [48l[54H56]. 

For a comprehensive discussion and literature survey of 
no-hair theorems, the reader is referred to the reviews of 
Bckenstein [57] and Chrusciel, Costa, and Heusler [58] . 
Although the literature is vast, there are two basic as- 
sumptions lying at the heart of most no-hair theorems. 
The first is that of stationarity, whose necessity has been 
demonstrated by Jacobson's "Miracle Hair Growth For- 
mula," a perturbative construction of a hairy BH with 
time-dependent scalar boundary conditions [5^ . 

The second assumption is the truncation of the scalar- 
tensor action to second order in the derivative expan- 
sion. At this level, the most general action [Eq. (pi) in 
the single-scalar case] is very simple and contains only 
three terms, whereas scalar-tensor gravity at the four- 
derivative level is too complicated to be studied in com- 
plete generality, and thus, attention is often restricted to 
particular models. 

One such model is quadratically modified gravity, 
whose action contains all possible terms quadratic in 
the Riemann tensor, coupled to a scalar. In this the- 
ory, BHs have been studied perturbatively, and solu- 
tions with scalar hair have been found [50^64] . More- 
over, in the special case of a scalar coupled to a topolog- 
ical invariant - namely, Einstein-dilaton-Gauss-Bonnet 
(EDGB) or Dynamical Chern-Simons (DCS) gravity - a 
no-hair theorem for neutron stars has been established 
[65j . The conclusion here is that spherically-symmetric 
neutron stars have vanishing scalar monopole moment, 
but higher-order scalar multipoles need not vanish. How- 
ever, the presence of derivatives higher than second order 
in the field equations severely complicates the implemen- 
tation of numerical simulations. 

Although four-derivative actions generically lead to 
field equations with more than two derivatives, there are 
some noteworthy exceptions. One is the Einstein-Skyrme 
(ES) system, a nonlinear sigma model with target-space 
5'[/(2), containing a term in the action quartic in scalar 
derivatives, and admitting linearly-stable BH solutions 
with scalar hair which have been described numerically 
in [55H55] . Another model with second-order field equa- 
tions is the galileon [69] , which is related to both higher- 
dimensional Lovelock gravity |70j and massive gravity 
[7T1 [7^ . It satisfies Solar System constraints by means 
of the Vainshtein mechanism |73j , and a no-galileon-hair 



theorem for spherically-symmetric BHs has been recently 
estabfished (74j . 



B. A generalized no- hair theorem 

The "classical" no-hair theorems described in the pre- 
ceeding paragraphs, which are statements about station- 
ary vacuum space-times, have been extended to the con- 
text of a BH binary system. Employing the "general- 
ized EIH" formalism developed by Eardley [75] , Will and 
Zaglauer [T7] have shown that the leading-order post- 
Newtonian (PN) dynamics of a BH binary in Brans-Dicke 
theory is indistinguishable from that in general relativity. 
Recently, this result has been extended to general scalar- 
tensor theories in the extreme mass-ratio limit [34' and 
it has been shown to hold up to 2.5PN order for generic 
mass ratio [76]. Thus, even a dynamical (vacuum) space- 
time with two interacting BHs does not have scalar hair 
in the PN limit. We can regard this conclusion as a "gen- 
eralized no-hair theorem." 

The generalized no-hair theorem relies on the assump- 
tion that the binary system is isolated, in the sense that 
cosmological and environmental effects (say, due to the 
galactic background surrounding the BHs) are neglected. 
More precisely, it is assumed that: (1) there is no matter 
in the system, (2) the scalar field has zero potential, (3) 
the scalar-tensor action is truncated to second order in 
the derivative expansion, and (4) the metric is asymp- 
totically fiat (in all conformal frames) and the scalar is 
asymptotically constant. 

The vacuum assumption (1) can be relaxed either by 
considering BHs in astrophysical environments or by con- 
sidering compact stars, which are affected by the well- 
known spontaneous scalarization phenomenon (cf. |77j 
for recent numerical studies). Other recent numerical 
studies created a scalar-field "bubble" around the bi- 
nary by using a nonvanishing potential, i.e. relaxing as- 
sumption (2). They found that the scalar-field bubble is 
rapidly accreted by the BHs, modifying the binary dy- 
namics [781. As for assumption (3), compact binary dy- 
namics in quadratically modified gravity has been studied 
analytically in a perturbative framework, where one takes 
the point of view that the model should be considered 
as an effective low-energy theory [55] [721 IHD] ■ Whether 
these theories are well posed for numerical evolutions is 
currently a matter of debate. We will not consider this 
problem in the present paper, but it is an interesting 
topic for future research. 

Relaxing either assumption (3) or assumption (4) in- 
troduces a new length (or time) scale in the BH binary 
dynamics. In the case of assumption (3) this scale is de- 
termined by the Compton wavelength of a heavy particle, 
whose square enters into coefficients of four-derivative 
terms in the action. In the case of assumption (4), the 
new scale is determined by cosmological and/or galactic 
effects. A priori, it is not obvious which of these effects is 
dominant, and thus it is worthwhile to explore both pos- 



sibilities. To our knowledge, the relaxation of assumption 
(4) has not been explored in the literature, and it is the 
main focus of our paper. 

Asymptotic flatness of the metric is only an approxima- 
tion to the dynamics of an astrophysical binary. Observa- 
tions show that the Universe is expanding on timescales 
which are very large, but nevertheless finite with re- 
spect to astrophysical BH binary evolution. As shown in 
[5^ [ST] , imposing time- varying boundary conditions en- 
dows the BHs in a binary with scalar charge, and there- 
fore the binary can emit dipole scalar radiation. Fur- 
thermore, many cosmological models consider the exis- 
tence of background scalars which can be anchored on 
matter [821 - I9T] . In this case, one can for instance con- 
ceive of a BH binary evolving in the background of a 
nearly-static but nonuniform scalar field anchored on the 
galactic matter. The characteristic lengthscale of such a 
scalar-field profile would be much larger than the binary 
separation, and therefore it would have the same effect on 
the dynamical evolution of the binary system as the en- 
forcement of boundary conditions which are not asymp- 
totically fiat. Finally, BH dynamics in the background 
of scalar fields could also be relevant to understanding 
accretion inside hypothetical supermassive boson stars, 
where huge scalar-field gradients are expected [92] . 

Scalar-field gradients can therefore allow us to cir- 
cumvent the generalized no-hair theorem, i.e., to have a 
spacetime which contains only BHs and still emits scalar 
radiation. Indeed, as we shall discuss below [cf. Eq. (28), 
Section IIIB2|, in the presence of a spatially- varying 



scalar-field profile ip{x), a nonrotating BH of mass M 
with world-line {t,x{t)) would have a scalar chargaH 



Q{t)=4M^^-V^{x{t)) 



(1) 



If the BH is accelerated, or if the scalar gradient is 
nonuniform, the scalar charge would evolve in time, yield- 
ing scalar radiation. 



C. Plan of the paper 

The main goal of this work is to explore the conse- 
quences of the presence of a scalar-field gradient, which 
corresponds to imposing nontrivial boundary conditions 
on the dynamics of a BH binary, and to verify numer- 
ically whether, as suggested by Eq. ([l]), such a setup 
can allow scalar radiation from a BH binary system in 
scalar-tensor theory. The plan of the paper is as fol- 
lows. In Section IH] we lay out our theoretical framework 



^ The scalar charge Q and mass M are Einstein-frame quantities, 
and geometrical units G = c = 1 are employed, where G is the 
Einstein- frame bare gravitational constant. The result quoted 
here also assumes that the BH motion relative to the scalar-field 
profile is "slow", in the sense that M{dx/di) ■ S/if <C 1. 



by introducing generic scalar-tensor theories and present- 
ing the relations that allow us to transform between the 
Jordan frame (where physical quantities should be com- 
puted) and the Einstein frame (where we will perform 
our calculations). In particular, we show how gravita- 
tional radiation in the Jordan frame can be computed 
from a knowledge of the Newman-Penrose scalars in the 
Einstein frame. In Sect ion [HT] we introduce analytical ap- 
proximations for scalar fields in the background of single 
and binary BH spacetimes, that will be used to validate 
(and provide insight into) our numerical simulations. In 
Section [TV] we present the details of our numerical imple- 
mentation. The results of our simulations are discussed 
in Section |V| In Section |VI| we summarize our findings 
and point out possible directions for future work. To 
improve readability, in the Appendices we collect techni- 
cal material that illustrates various important points of 
our analysis. Appendix [A| shows that a BH moving with 
constant velocity in a uniform scalar-field gradient does 
not emit scalar radiation. Appendix ^ (which is com- 



plementary to Section HI B 2 ) collects some lengthy for- 



mulas illustrating the structure of gravitational radiation 
from a BH binary moving in a scalar-field gradient. In 
Appendix[C]we provide explicit expressions for the evolu- 
tion equations used in our numerical code. In Appendix 
[D] we estimate the order of magnitude of the scalar-field 
gradients expected in scalar-field models of dark matter. 



II. THEORETICAL FRAMEWORK 

We focus on general single-scalar-tensor theories in 
vacuum with vanishing scalar potential, and at most two 
derivatives in the action. These theories are equivalent 
to Einstein's theory extended to include a minimally cou- 
pled scalar field with vanishing potential. This statement 
(which will be clarified below) has a nontrivial conse- 
quence: the addition of minimally coupled scalars to Ein- 
stein's gravity allows one to study a multitude of scalar- 
tensor theories at once. For this reason our simple frame- 
work offers an opportunity to take a glimpse at a rather 
broad spectrum of physics beyond Einstein's theory. 

Our starting point is the action of a general scalar- 
tensor theory for a single scalar field (p, written as 



S 



d^x 



-9 



16ttG 



(F(0)i? - 87TGZi(^)g^''d^(f>d,P - 



u{4>)) , 

. (2) 
where R is the Ricci scalar associated to the metric g^^^ 
and F {(/)), Z {(p) and C/(0) are arbitrary functions (see 
e.g. [H and references therein). This form of the ac- 
tion corresponds to the choice of the so-called "Jordan 
frame" , where the scalar field is nonminimally coupled 
with gravity and all other matter fields obey the equiva- 
lence principle. The dynamics of matter fields would be 
described by an additional term Smatt on the right-hand 
side, that we set equal to zero because we are interested 
in BHs in vacuum. 



Our numerical evolutions are more easily performed 
in the so-called Einstein-frame representation, that is re- 
lated to the Jordan-frame representation by a conformal 
rescaling of the metric. In the Einstein frame, the scalar 
field is minimally coupled with the metric tensor and 
it affects the scalar-matter coupling in the matter ac- 
tion S'matt- Working in the Einstein frame is convenient 
because we focus on pure BH spacetimes, i.e., we set 
•Smatt = 0. Since we are interested in the effect of bound- 
ary conditions, we shall assume for simplicity (as in [8T] . 
and at variance with [75]) that the effect of the scalar- 
field potential is negligible: C/((/)) = 0. 



A. From Jordan to Einstein and back 

With U{(j)) — 0, the explicit transformations that re- 
cast the previous action in the Einstein frame are [20] 



5^. = FWa^ 


IV 1 




(3) 


^i(t>) = Jd^ 


'iF'{cf)f SnGZicj))' 
[2F{cl>)^ ' i^(0) J 


1/2 


(4) 


A(^)^F-i/2 


0). 




(5) 



The Einstein-frame action is then 



5 = Y^ / [i?"^ - a'^Ed.vd^v] V^d^x , (6) 



and it leads to the following equations of motion: 

1 



G" 



^^ - df,(pd^(p - ^g'^^g''j^dpipdaV , 



n^(^ = 0. 



(7) 
(8) 



where the label E denotes quantities built out of the 
Einstein- frame metric g^^. This will be the starting 
point of our analysis. It is important to stress that, even 
though we are formally investigating a minimally coupled 
theory, our results are in principle applicable to a wide 
range of scalar-tensor theories. By working in the Ein- 
stein frame we can focus on quantities that depend on the 
intrinsic properties of the binary, rather than quantities 
that would be measured by gravitational- wave detectors. 
The latter may be obtained for any specific theory using 
the transformation from the Einstein frame to the Jordan 
frame. 



In the Einstein frame this corresponds to (7^^ 



and from Eqs. P^-p^ one gets 
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(1) 



(10) 



(11) 



In general, gravitational waves in scalar-tensor theories 
have 3 degrees of freedom [321 IM] ■ In the Einstein frame 
they correspond to the two transverse-traceless compo- 
nents of the metric perturbation, plus the scalar field. 
The calculation of these quantities does not present dif- 
ficulties. The physical degrees of freedom, which should 
be computed in the Jordan frame, can be read off from 
Eqs. (lO)-(ll). An alternative procedure is presented in 
Ref. |77| . which shows the transformation of the corre- 
sponding Newman-Penrose quantities, with similar re- 
sults. In this work we will present results for the scalar 
field (^ and for the curvature scalar ^f" in the Einstein 
frame, which is directly related to the two Einstcin-frame 
polarization states h^ and h^ via 



y^f ^h%- ih^ , 



(12) 



where dots denote time derivatives. Quantities in the 
Jordan frame can be found using Eqs. (lO)-(ll), once 



one specifies the underlying theory. From here onwards, 
having established the relation between the Einstein and 
Jordan frames, we shall work exclusively in the Einstein 
frame, dropping the label "E" from all quantities. Unless 
specified otherwise, we will use geometrical units and set 
G = c = 1. Note however that here G is a bare grav- 
itational constant, and it is different from the quantity 
measured by a Cavendish experiment. 



III. SCALAR FIELDS IN SINGLE AND 

BINARY BLACK-HOLE BACKGROUNDS: 

ANALYTICAL APPROXIMATIONS 

In this Section we introduce approximate analytical 
solutions that describe single and binary BHs in a scalar 
field gradient. These solutions will be useful below, either 
as code checks or for the interpretation of our numerical 
results. 



A. Single black holes: 
linearized analytical solutions 



B. Gravitational w^aves in the Einstein and Jordan 
frames 



In the Jordan frame, we describe perturbations of the 
metric and of the scalar field as follows: 



9^. = gl^J + V , </> = <t>^'^ + ^^^^ 



(9) 



Let us assume that the scalar-field gradient is of such 
low amplitude that the scalar can be treated as a per- 
turbative effect on the spacetime metric. Under this as- 
sumption we can neglect terms quadratic in the scalar 
field, and therefore the field equations reduce to 



Rfiv — , 



(13) 
(14) 



We will drop this perturbative approximation in Section 
IV where the field equations ^ and (Is]) will be solved 
numerically. 

Equation ( [l3| is of course identical to Einstein's equa- 
tions in vacuum. We will consider the Schwarzschild (and 
later the Kerr) metrics as background BH solutions, and 
we will solve the Klein-Gordon equation (14) on these 
backgrounds. For the reasons explained in the intro- 
duction we are interested in numerical evolutions in a 
scalar-field gradient. Therefore we will consider back- 
ground scalar-field solutions generated by distant, fixed, 
infinite homogeneous planes with constant surface scalar- 
charge density a. To our knowledge, Press [US] was the 
first to study a closely related problem (his setup differs 
from ours in that he considered a spherical shell of scalar 
charge as the source of the scalar field) . Here we recast 
some of his results in a form suitable for comparison with 
our numerical setup. 



1. Spherically symmetric black-hole background 

Let us first consider the Schwarzschild solution in 
isotropic coordinates: 



not emit scalar waves. Indeed, as we will show analyti- 
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related to the areal radius r by r = f (l 



(15) 



z^ is the isotropic radius, which is 



M-) 



The scalar-field equation ( 14 ) on this background ad- 



mits a simple axisymmetric solution sourced by distant 
planes of constant scalar-charge density, i.e.: 



<Poxt 



27ro'(r — M) cos ( 
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cos 9 
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4f2 J ' 



(16) 



where z = f cos 9 is the direction orthogonal to the 
charged plane. If no BH is present, Eq. (161 reduces 
to the field tpcxt = 2'Kaz generated by homogeneously 
charged infinite planes, with constant gradient dzfcxt — 
2TTa. For large \z\ the constant-gradient behavior applies 
also to the case where the background is a Schwarzschild 
BH. 

We shall take as initial condition ipini = 2'kijz and test 
our numerical framework by checking that, after a tran- 
sient, the scalar-field profile settles to the analytical so- 
lution </j = (^cxt , up to corrections of second order in the 
scalar field. 

In Appendix |A] we derive a solution describing a BH 
moving with small constant velocity in a direction or- 
thogonal to the charged planes, and show that it does 



cally in Section HI B and numerically in Section IV the 
BH must have nonvanishing acceleration in order to gen- 
erate scalar radiation. 



2. Rotating black-hole background 

In the case of a rotating BH with a rotation axis that 
is not orthogonal to the charged planes, axisymmetry is 
lost; however, a simple solution can still be found |95) . 
Let us choose our coordinates so that the BH spins along 
the z-axis, at an angle 7 with respect to the direction 
orthogonal to the charge-carrying planes. Then it can be 
shown that a solution is 



ifcxt — 2TTa{r — M) (cos 7 cos 9 + sin 7 sin 9 cos 0/q) 



27rcr(r - M) 
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J a 



r- M 



7^ 



r[2 - la/S] p,a,s f r-M 



(17) 



, (18) 



r denotes the F- function and P!^{z) is an associated 
Legendre function of the first kind. At large distances 
/a — >■ 1, and it can be seen that the charged plane gener- 
ates a uniform gradient in the xz-plane. Figure [l] shows 
contour plots of (^cxt in the y = plane, for different 
values of a and of the angle 7. Notice how the field lines 
are distorted and frame-dragged close to the horizon. 



B. Black-hole binaries: analytical approximation 
for quasi-circular inspirals 

Let us now turn to the more complicated case of a 
quasi-circular BH binary evolving in an external scalar 
field with constant gradient. We begin our discussion by 
addressing the delicate problem of specifying boundary 
conditions for this system, which will be used to perform 
the numerical simulations discussed in Sec. IIVI 



1. Boundary conditions 

In our numerical setup the boundary conditions are im- 
posed by fictitious faraway charges, which are not mod- 
eled in our numerical simulations. These charges are as- 
sumed to lie outside the numerical grid, and mimic an 
external profile due (say) to a galactic scalar-field back- 
ground. This "large-scale field" acts as a sort of reser- 
voir: when the local field begins falling into the BH, there 
is an ingoing flux which restores the scalar-field gradi- 
ent. The presence of this field outside the boundaries 
of our numerical simulations is implemented through the 
boundary conditions. At large distances we want to al- 
low for outgoing waves while imposing the existence of 




FIG. 1. Contour plots of the field i^oxt in the vicinity of a rotating BH, as given by Eq. ( 171. Top: The infinite charged plane 
is at an angle 7 = 0, and the BH has dimensionless spin a — (left) and a/M = 0.99 (right). The value of ipcxt/{2TTa) is 
shown along selected contour lines; the two panels only differ because of the different size of the horizon. Bottom: A BH with 
a/M = 0.99 is immersed in a field gradient at angles 7r/4 (left) and it/2 (right). All contour plots refer to the plane y = 0. 
Selected contour lines correspond to the same values as the top panels. 



the scalar-field gradient, so we require the field to behave 
in the following way: 



the homogeneous equation Dtp = describing outgoing 
(approximately spherical) scalar waves. We thus get 



</5 = </5oxt + 



$(i 



(19) 



where r is the areal radius, the constant gradient corre- 
sponds to the external field (/?oxt = 2Tra{r — M) cos 6, and 
the second term on the right-hand side is the solution of 
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(rifie^t) 



(20) 



Since the boundary conditions are defined at large dis- 
tances, ifcxt — '2nar cos 9, and we can write 



dr 



d 

{rip) + — {np) — Airar cos 9 . (21) 



2. Multipole expansion of the scalar field 

The angular dependence of the scalar field can be de- 
scribed through a multipole expansion of the form 

1 
2' 



$(t-' 



= M + n'V, 



-n'n^Qij 



(22) 

where $ is the function appearing on the right-hand side 
of Eq. (191, dots denote derivatives with respect to the 
retarded null coordinate u = t — r and n — x/r = x/\x\ is 
the radial unit vector, which depends only on the angles 
{9,(j))- The calligraphic symbols Ai, Vi and Qij denote 
the monopolar, dipolar and quadrupolar components of 
$, respectively. 

The full relativistic scalar equation to be solved is 
DiyS — 0. In order to obtain a solution which describes 
the physics that we are interested in - namely, a BH 
binary in a scalar gradient - it is essential to impose cor- 
rect boundary conditions, both at null infinity, and in the 
vicinity of the worldlines of the singularities of the two 
BHs. The boundary conditions at null infinity have been 
discussed above; the near-worldline boundary conditions 
that must be imposed in order to find an approximate 
solution are more complicated. 

In the special case of a comparable-mass BH binary 
system with small size-to-separation ratio (or, equiva- 
lently, small orbital velocity) the problem of imposing 
correct near-worldline boundary conditions is substan- 
tially simplified when one employs a "point-particle" ef- 
fective field theory, in which length scales smaller than 
the Schwarzschild radii are integrated out. In this effec- 
tive field theory the matter action has the form 



CPP _ \^ 



dsA^A , 



(23) 



where A is an index that runs over the bodies {A — 1,2 
for a binary system), Ta is the worldline of body A, dsA 
is the proper differential arclength along Ta, and Ca is 
the "effective point-particle Lagrangian" of body A. For 
a structureless particle of mass uia, one has Ca = —ruA- 
The matter action ( 23 1 gives rise to sources in the field 



equations of the effective point-particle theory, and these 
sources automatically enforce the correct boundary con- 
ditions at the worldlines Ta for both g^^, and (p. For 
instance, to leading nonrelativistic order, the scalar-field 
equation has the explicit form 



Dfp = 4np^ 



= 47:J2QAS^^\x~ZA{t)), 



(24) 



where D/ is the D'Alembert operator in flat space, 
(tjZAit)) is an explicit parametrization of the worldline 
r^, and Qa are the scalar charges of the BHs. 

Moreover, to leading nonrelativistic order, one finds 
that the multipole moments entering into Eq. ( 22 ) are 
given by 



(25) 



M = / d-^xp^ = Qi+Q2 



V 



d XX p^ 



QlZl + Q2Z2 ; 



(26) 



and so on. 

In general, calculating the scalar charges Qa is a dif- 
ficult problem. A simplification takes place if we assume 
that the BH masses Ma have the same order of magni- 
tude M, and that 



cr< 



a 
AP 



1 

Maj 



(27) 



where a is the typical orbital separation, and v is the 
typical orbital velocity. Henceforth, terms of order 
(Afwcr)^ will be dropped. Then the BH scalar charges 
Qa rnay be found by Jacobson's formula [59 , which for 
Schwarzschild BHs yields 

'dip{t,ZA{t)) 



i(.t) 



AMI 



dt 



VA{t)-Vip{t,ZA{t)) 



^SiraMlvAit)- z, (28) 

where VA{t) = ZA{t) is the velocity of body A, and in the 
second line the full scalar field (p{t, x) has been replaced 
by the zeroth-order field (^cxt = 27rcrz. 

Let us specialize to quasi-circular orbits in the yz- 
plane, so that the trajectories take the simple form 

Zrc\{t) = Zi{t) - Z2{t) 

^a{t)[y COB x{t) + zsinx{t)], (29) 
Mzcmit) = M^z^{t) + M2Z2{t) , (30) 

where a{t) is the orbital radius, x(0 — J uj{t)dt is the 
orbital phase, uj(t) is the angular frequency of the orbit, 
and zcuit) is the center of mass of the binary system. 
The quantities d, x — uj, and ocm = zqm are all small (of 
order 1/c^), and vanish in the absence of radiation reac- 
tion. An explicit expression for their leading-order time 
evolution may be found by solving the 2.5PN equations 
of motion (given in Section 9 of I96J), while dropping 
conservative corrections. Carrying out this calculation 
yields 

256ty/''_ 
256t\^^/* 



a{t) = a(") (1 - 



,(0) 



64 A 

5^y 



(31) 



m = x^"^ 1 1 



STqj 



^(0) 



96A 

5^y 



(32) 



where a^^' and x sue the (constant) radius and angular 
frequency of the zeroth-order orbit, respectively, and 

\-\ -1 



A=l 



MiM2{Mi+ M2) 
[a(0)]4 



(33) 



is the time scale over which the quadrupole tensor radi- 
ation shrinks the orbit. 

With an exphcit description of the orbit in hand, the 
scalar charges Q1.2 and multipole moments ( 25 )-( 26 ) may 
be calculated (see Appendix Ib] for the explicit expres- 
sions). In this way, we find that: 

1) monopole radiation is emitted at the orbital fre- 
quency: 



M 



Sttct 



M{Ml+Ml){vcM-z) 



(34) 



+ MiM2{Mi - M2)(asinx + axcosx) 



and it vanishes in the equal-mass limit; 

2) dipole radiation is emitted at twice the orbital fre- 
quency, and more precisely 



Therefore the m = contribution should oscillate with 
frequency $7, the to = 1 contribution with frequency 2Q, 
the m — 2 contribution with frequencies Sfi and fi, and 
so on. As we will show below, this behavior is consistent 
with our numerical simulations. 



IV. NUMERICAL IMPLEMENTATION 

Our numerical implementation of scalar-tensor theory 
closely parallels [57 . but borrowing notation and con- 
ventions from [US]. The physical system studied in [51] 
was quite different, since that paper considered higher- 
dimensional Einstein gravity in vacuum. However the 
equations can be cast as a system involving a scalar field 
coupled to gravity via dimensional reduction, so they are 
formally similar to the system considered here, as we 
show below. 



v^n 



CM 



2?rcl, DC + Acl, c 



(35) 



where the "CM" term is emitted at the orbital fre- 
quency, the "DC" component is nonoscillatory, and 
the I?ioi, osc component oscillates at twice the or- 



bital frequency: cf. Eq. ( B9 1 



The physical problem addressed here differs from the 
situation investigated in [STj. That study considered 
time-dependent scalar boundary conditions, rather than 
a gradient, and it found that monopole radiation is ab- 
sent, while dipole radiation vanishes in the equal-mass 
limit. One may expect dipole scalar radiation to be emit- 
ted at the orbital frequency, rather than twice the orbital 
frequency. The reason why this expectation is erroneous 
in our case is that we have a background field with a gra- 
dient directed along the orbital plane, which combines 
with the oscillatory component sourced by the orbital 
motion. 

A simple toy model can provide us with a comple- 
mentary and perhaps more intuitive way to justify the 
expectation that dipole radiation must be emitted at 
twice the orbital frequency. Let us consider a rotating 
source with frequency fi on a scalar-field background 
Vext = 27r(T2 = 27rcrr sin 9 sin (f) [in our "rotated" polar 
coordinates, see Eqs. (74) below]. The source will pro- 



duce a modulation in the background field of the form 



<f^^,^t[l + f{cl)-nt)]. (36) 

Expanding in circular harmonics, f{4> — Qt) — 
E™/me™(*-"*)and 



(p = 27Tar sin 6 sin 0(1 + > /„ 



^i?n(0— Sli) 



), (37) 



which implies that the multipolar components of the field 
will have the following dependence: 

Vim ~ (e-'(™+i'"* + e-i(™-i)S") + constant . (38) 



A. 3-1-1 decomposition 

As a preliminary step for our numerical implementa- 
tion, we perform a 3 -t- 1 decomposition of the spacetime 
(see [55] and references therein). Let us consider a shc- 
ing of the spacetime in a set of three-dimensional surfaces 
S. Introducing the normal n^ to the surface S and the 
projector 



Tm^ — 9l^tl^ I Ti^Tlu , 



(39) 



we write the four-dimensional metric in the form (/i, v = 
0,...,3;i,j = 1,2,3): 

ds = g^j,^dx^dx" 

= -a^dt^ -t- -i,j{dx' + (3'dt){dx^ + p'dt) , (40) 

where a, /3* are the lapse and the shift, respectively, and 

dt = an + l3. (41) 

We shall denote by Di the covariant derivative on E, 
i.e., the covariant derivative with respect to the three- 
dimensional metric 7^-. 

The Lie derivative with respect to n^, then, is £„ = 
[dt — Hfi)l a. Defining the extrinsic curvature 



-f^ij — n'—nlij 



(42) 



we have the following evolution equations for the (three- 
dimensional) metric: 



{dt - Cp)^i.j = -2aKi 



Furthermore, we define the scalar curvature K^ as 



K,, 



-Cn'P 



(43) 



(44) 



so that 



{dt - Cp)ip 



-2aK,, 



(45) 



Comparing with the definitions of the variables Qi and 
n in I97j, we find that 



Since n''df,(p = ~2K^, ^''"d^ip = D^ip, g^^ 
n^n^, defining, as on page 87 of [99] 



7^1^ 



f^i^ ) 



n 



Dup _ -/if.df'ip 



we get 






Therefore, 



K^ ^ -^VStK^H . 



(46) 



(47) 



BttG p = n^n''df,ipdyip - -gf,uda'pd°'ip 



2Kl + -D,^D\, 






(55) 
(56) 



(57) 

(58) 

(59) 
(60) 



The constraint equations (2.15) and (2.16) of (HZI read 

n2 I n2 



f ^ 



D.^D'if, 



where we have used the fact that 7°*^ = (see |99l), thus 
j^'^Di.ipD^'f = j'W.ipDjip. Therefore, Eqs. (2.4^.6) and 
(2.4.9) coincide with our constraint equations ( [48| and 
(49 1. Furthermore, we can compute the quantity (cf. 



D,K\ 



D,K = 



no, 

/ 



= 2K^D,ip . 



(48) 
(49) 



page 89 of [Ml) 



We have 



5^" = 7''"7"'^T„,3 . 



(61) 



The evolution equation (2.17) of [97] can be written as 
= -D^Dja + a (^^''Rij + KK,j - D,(pDjip - 2K,kK 



SttGS'"' = D^^D^if - -j'^'D'^D.ip + 2Y''Kl , (62) 



fe 
(50) 



and the trace of this equation (since 7^ = 3) yields 
1 



Note that since ■y'^iJ^jR^u — DifDjip, this expression 
coincides with Eq. (2.23) of ^ 
Eq. ( 50 ) we have 



i7rGS^--D'^D,ip + 6K'^. 



Taking the trace of Then, 87tG{S - p) = 4K^ ~ D.ipD'ip, and 

AttG [{S - p)7,, - 2S,,] = -D,ipDjip 



JV^' 



(63) 



(64) 



dtK-P^diK+D'Dia-a {^'^^ R + K^ - D\D,ip^ , (51) therefore Eq. (2.5.6) of [99 coincides with our Eq. ^ 



and using Eq. ( 48 1 we find 



{dt - Cp)K = ~D'D,a + aK.jK'^ + AaK^ , (52) 

that should be compared to Eq. (2.19) of [HZj. The evo- 
lution equation (2.18) of [97] can be written as 



1 



{dt-Cf3)K^^-- 



1 



;£„n 



= K^K D'ifD.a ~ -D.D'ip . (53) 

2a 2 

All terms of this expression appear, with the same coeffi- 
cients, in Eq. (2.37) of [SSj. The additional terms in that 
equation which are not present here are due to the more 
complicated dynamics of the scalar field arising from di- 
mensional reduction. 

It can be useful to write the equations also in terms of 
the stress-energy tensor, which enables us to compare the 
scalar- field terms with those in 99J. From (It]) we have 

8ttGT^„ = df^ipd^ip - -gu.^da<pd°'ip . (54) 



B. Baumgarte-Shapiro-Shibata-Nakamura 
formalism 



Our evolution equations use the Baumgarte-Shapiro- 
Shibata-Nakamura (BSSN) formalism, in which the dy- 
namical variables are {x,7ij, ^y, -f^, F*}, defined as fol- 
lows: 

lij = X^^lij (with 7'^ = xf^) , 
X=(det7y)-i/3, 



A J - X [ K,j - -7„X 



F'==y^Ff,-xf^ + ;^7'^aa- 



(65) 



Alternative notations replace our variable x by a vari- 
able "0 defined as x~^ — V"^- Here we will use x ^^ our 
dynamical quantity. 

The Einstein equations in the BSSN formulation in 
the presence of a scalar field, which appears through the 



10 



quantities p,j^, Sij, are implemented in an extended ver- 
sion of the Cactus [lOOj based Lean code [IQll 1102 
in the form given in Eqs. (Cl)-(CIO) of Appendix IC 
Mesh refinement for our simulations is provided by Car- 
pet 103 ■ horizon diagnostics by AHFinderDirect 



[104| I105J and BH binary initial data satisfying Eq. ( 13 ) 
by a spectral solver |106j provided through Cactus as 
the TwoPuNCTURES thorn. 



C. Initial data 

Our initial data sets consist of either a single BH or 
a BH binary evolving in a background scalar field with 
a nonvanishing gradient. The background scalar field is 
generated by distant sources, which are kept fixed in our 
time evolutions. As a simple way to enforce this scenario, 
imagine that the background scalar field is generated by 
infinite homogeneous charged planes with surface den- 
sity a. As we saw in Section |III[ if a BH is present in 
the spacetime and the scalar field is small enough to be 
treated in the linear approximation the metric is unaf- 
fected, but the equilibrium solution for the scalar field 
changes. This setup is quite similar in spirit to the one 
adopted by Palenzuela and collaborators [l OTl 1108] , the 
main difference being that they dealt with electromag- 
netic fields, and that their external magnetic field is gen- 
erated by a current loop far away from the system. 

Under these assumptions, the solution (^cxt we found 
in Eq. ( 16 1 should be a good approximation to the initial 



data describing a single, nonrotating BH in a scalar-field 
gradient. Except for regions close to the BH singularity, 
the term M^/4f^ in parentheses in Eq. (16) can further- 
more be neglected. Therefore we initialize the scalar field 
using the simplified expression ipcxt = li^az. Our numer- 
ical simulations of single BHs show indeed that the initial 
data ( [l6| and its approximate version 2t:gz yield (after 
a brief transient, which we exclude from our analysis be- 
low) virtually identical evolutions of the scalar field and 
of the spacetime metric. 



D. Upper bounds on the field gradient from the 
threshold of black-hole formation 



Our setup consists of a constant scalar-field gradient 
at large distances. Because the energy density in any 
classical field theory is proportional to the square of the 
field gradient, one expects a roughly constant energy 
density p ^ a^ . The total mass M in a region of lin- 
ear dimension R then scales like ~ pR? , and therefore 
MjR ^ pR^ ~ a'^R^. Therefore we expect that the ini- 
tial data will contain a horizon for a > R~^. This con- 
dition imposes a nontrivial constraint on the size of our 
numerical grid. Here we provide a more formal argument 
supporting this conclusion. 

We focus on conformally flat backgrounds with A^ i = 
0, H = 0. Then the momentum constraint is identically 



satisfied, while the Hamiltonian constraint yields 

= A-0 + -i/jrf^dtipdjip , 



(66) 



where rj^^ is the Minkowski metric. A scalar field with 
constant gradient is such that ip — lixaz^ and therefore 
we get the equation 



A^ 



7rV2 



V'- 



(67) 



Imposing regularity at r = 0, the solution of this equation 
is 



i, 



A sin [i^ar l\Ji) 



(68) 



An apparent horizon exists for (?/'*r'^)' = 0. This condi- 
tion is equivalent to finding the roots of 



a; cot a; — 1/2 , 



ith X = 



72 



(69) 



The smallest root of this equation is at x ~ 1.16556, i.e. 

, ^/2 0.52469 



cr~ 1.16556- 



(70) 



in good agreement with the previous back-of-the- 
envelope estimate. 

Our boundary conditions are enforced at distance 
r jM = 160 and 384, respectively, for single and binary 
BH simulations. This provides a formal upper bound of 
Ma w 5 X 10"'^ on the magnitude of the gradients that 
we can simulate. However, as we will see below, even 
smaller gradients Ma « 10^* can generate exponentially 
growing instabilities (indicating collapse) in our numeri- 
cal evolutions. 



V. NUMERICAL RESULTS 

In this Section we first discuss our results for single 
BH evolutions, verifying that for small values of the gra- 
dient Ma they are in good agreement with the analytical 
predictions of Section III A Then we study gravitational 
and scalar radiation from BH binaries in a scalar-field 
gradient. 



A. Single black-hole evolutions 

In order to compare our numerical results with the 
analytical predictions of Section |III[ it is important to 
minimize coordinate effects. Let us denote by f the ra- 
dial coordinate used in our numerical simulation, which 
coincides with the isotropic coordinate at i = 0, i.e. 
f{t = 0) = f. This need not be true at later times, since 
the gauge can dynamically change during the evolution. 
In order to monitor the scalar field as a function of time 



11 



0.010 
0.008 
„0.006 
"^ 0.004 
0.002 
0.000, 



- 


' 1 


error (%) 


IB 
6 
4 
2 



ri 1 'u 


10 29 3Q 40 50 


- 








r/M 




















- 






























- 






200 400 600 800 1000 1200 

t/M 




200 300 

t/M 



500 



FIG. 2. Real part of the scalar dipole mode </5io (the imagi- 
nary part vanishes) for Ma — 10"^ and extraction radii (from 
top to bottom) f/M = 50, 40, 30, 20, 15, 10 and 5, compared 
to the predictions of Eq. (76 1. Solid lines refer to the numer- 



ical evolution; dashed lines refer to the analytical solution 
evaluated at time-dependent areal radii r, which are com- 
puted dynamically during the evolution. The inset shows the 
percentage discrepancy between the numerical and analytical 
prediction as a function of extraction radius. 



and check that it eventually settles to a configuration 



close to the analytical solution (/Sext of Eq. (16), we need 



to perform integrations over spheres at given values of f . 
Furthermore, we compute the areal radii r at these loca- 
tions by performing spherical integrations of the metric 
components, as discussed in Ref. 



The spherical-harmonic cxpansiorjj of the scalar field 
reads 

^{t, r, e,q^)^J2 '^'™(^' r)Yi,n{e, 0) , (72) 



Itn 



where 



^i™(i,r)= dncpit,r,e,^)Y;^{9,^)^ 



(73) 



Since the binary moves in the yz-plane, we use noncanon- 
ical rotated coordinates 



X — r cos a , 
y = r sin 9 cos ( 
z = r sin 9 sin o 



(74) 



^ This multipolc expansion can be simply related to that intro- 
duced in Eq. |[22l. For example, we have: 



V-. 



ipiox - 



Vil ~ Vi-l - 



.ipil -I- '/'I -1 . 

I — z 

V2 



(71) 



l/r 



where the subscript l/r denotes that only the l/r term should 
be kept. 



FIG. 3. Real part of the scalar dipole mode ipio (rescaled 
by Ala) at the largest extraction radius f/AI = 50 for 
Ma = 10 ■ 
of Eq. 



* and Ma = 10 '^j compared to the predictions 

( 76 1 . Solid lines refer to the numerical evolution; 



dashed lines refer to the analytical solution evaluated at time- 
dependent areal radii r, which are computed dynamically dur- 
ing the evolution. The evolution does not settle to the ana- 
lytical solution for Ma = 10"**: there is an exponentially 
growing mode. This also shows up as an exponential growth 
of the subleading multipoles, as can be seen in Fig. Ill 



For small scalar-field gradients, the back-reaction on 
the metric is very small and we expect to recover the 
stationary solution (yScxt found in Section III A[ i.e., we 
expect that after an initial transient 



fit 



J, I , 1/ , V 



tPcxt [r, 9) = ipo + 27rcr(r - M) cos ( 



(75) 
where ipo is a constant. Because |Yio| = •y/3/(47r) cos6', 
the dominant nonvanishing component ipim at late times 
(up to quadratic corrections in the field) should be given 

by 



(76) 



<PlO = \/y27rcr(r-M). 



In this expression the areal radius r is effectively time- 
dependent, and it is computed dynamically during the 
evolution as described above. 

The numerical and analytical predictions (solid and 
dashed lines, respectively) are compared in Figs. [5] and 
|3] Let us focus first on Fig. [2] which refers to Ma = 
10^^. In this case the numerically extracted dipole mode 
asymptotes quickly to the analytical prediction. The in- 
set shows the percentage difference between the analyt- 
ical and numerical value of the scalar field at late times 
as a function of the extraction radius: the agreement is 
remarkable at large radii but it gets progressively worse 
as we get closer to the BH, most likely because of gauge 
effects (we can exclude that the deviations are due to non- 
linear effects, because these would scale with a^ , whereas 
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FIG. 4. Absolute value of the real part of the scalar multipoles |Re(iy9;o)| evaluated at the largest extraction radius f — 50M for 
different values of I and two values of the scalar-field gradient, Ma = 10~^ and A4a = 10~* (left and right panel, respectively). 



the disagreement seems to be independent of a) . To sum- 
marize, the analytical and numerical predictions agree 
within a few percent for gradients Ma < 10~^. 

The situation changes significantly for larger values of 
Ma. In Fig. |3] we overplot the numerical and analyti- 
cal components for Ma = 10^'' and for a larger gradient, 
Ma = 10"''. As expected, when rescaled by their respec- 
tive gradients the dipolar components are essentially the 
same, and they also converge to the value predicted by 
the analytical solution. However for Ma — 10~^ this con- 
vergence can only be seen at intermediate times, whereas 
at late times the mode develops an instability: the figure 
shows that the roughly constant "baseline" value of the 
field seems to be superimposed to an exponentially grow- 
ing oscillation that develops on timescales t ^ lO^M, 
both when we evaluate the field numerically (continuous 
black line) and when we use the areal radii to compute 
the analytical solution (dashed black line). 

The existence of this instability is confirmed by Fig. |4j 
where we look at higher multipoles of the field on a 
logarithmic scale for Ma = 10"^ (left panel) and for 
Ma = 10"'' (right panel). The right panel of this plot 
shows that for the larger gradient Ma = 10"* an expo- 
nentially growing instability with growth time ~ 10^ M 
is present also in the subdominant I = 0, I = 2 and 
I = 3 multipoles. The left panel illustrates that when 
Ma = 10"^ the instability (if it is present at all) de- 
velops on much longer timescales t > lO'^M, so it does 
not affect our numerical simulations. Notice that early- 
intermediate time results are physically consistent: for 
both values of the gradient the scalar-field distribution is 
dominated by the dipolar component, and it is in good 
agreement with analytical predictions. 

In summary, single-BH simulations in a scalar-field 
gradient show that our numerical evolutions are stable 
and reliable as long as the gradient is not too large. This 



conclusion is compatible with the arguments presented 
in Section flV Dl above. 



B. Binary black-hole evolutions: scalar and 
gravitational radiation 

In this Section we discuss our numerical evolutions of 
BH binaries in a scalar-field gradient. In this initial study 
we focus on evolutions of nonspinning, unequal-mass bi- 
naries with a gradient Ma = 2 x 10"^ and mass ratio 
q = 3, because unequal-mass binaries display two inter- 
esting features which would be absent by construction in 
the equal-mass case: center-of-mass recoil [see Eqs. (Bl)- 
(B3)] and monopole scalar radiation [see Eq. ( |B4| ]. 

Gravitational waveforms, as characterized by the 
Newman-Penrose scalar ^>4, are shown in Fig. [5] For 
such low values of the scalar-field gradient, the impact 
on gravitational radiation emission is hardly noticeable. 

The emission of scalar radiation is much more inter- 
esting. The scalar field acquires a nontrivial profile due 
to the dynamics of the orbiting BH binary. The scalar 
radiation has a crucial dependence on the binary setup, 
and more specifically on the angle between the orbital an- 
gular momentum of the binary and the direction of the 
scalar-field gradient. If this angle is zero, then effectively 
the individual BHs do not traverse any field gradient, and 
the scalar profile is expected to be trivial. Our numeri- 
cal results confirm this expectation: the output of these 
simulations is indistinguishable from vacuum evolutions 
in pure general relativity [1021 IllOj . 

On the other hand, the induced scalar radiation should 
be maximized when the orbital angular momentum is 
perpendicular to the field gradient, so we now focus on 
this case. Our results are summarized in Figs. [5] and [61 

Because the binary evolves on the background of a 
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FIG. 5. Numerical results for a BH binary inspiralling in a scalar field gradient, with the orbital angular momentum perpendic- 
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FIG. 6. Numerical results for a BH binary inspiralling in a scalar field gradient, with the orbital angular momentum perpen- 
dicular to the gradient. Left: dependence of the various components of the scalar radiation Ke{ipii)/{Ma) on the extraction 
radius (top to bottom: 112A'/ to 56M in equidistant steps). The dashed line corresponds instead to 10~'^lm{(pii)/{Ma) at 
the largest extraction radius. Right: time-derivative of the scalar field at the largest and smallest extraction radii, rescaled 
by radius and shifted in time. Notice how the two waveforms show a clean and typical merger pattern, and that they overlap 
showing that the field scales to good approximation as 1/f. 



dipolar scalar-field profile, this constant "background" 
value of the scalar shows up as a large imaginary com- 
ponenlP] of the I = \m\ = 1 scalar-field modes, which 
is apparent in the left panel of Fig. [6] (in fact, we had 



^ In both the singlc-BH solution | |l7| l with 7 = n/2 and in the 
numerical solution discussed here, the polar axis (in terms of 
which the polar angles, and then the harmonic decomposition, 
are defined) is orthogonal to the gradient. However the ordering 



to rescale the imaginary component by a factor 10^ in 
order to show this on this plot). 

The real part of ipn displays interesting dynamics (the 
imaginary component also has similar dynamics, but this 
is partially masked by the large background dipolar field. 



of the axes is different in the two cases. This explains why the 
imaginary part Im((/3ii) of the numerical solution corresponds to 
the real part Ke(ipii) of the analytical solution. 
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so the analysis of the real part turns out to be numeri- 
cally "cleaner.") At any extraction radius Re(</?ii) is ini- 
tially zero, as the binary is simply traversing a constant 
scalar field. As the binary evolves, we expect to see scalar 
radiation crossing the extraction surface and producing 
a nonvanishing scalar profile. This is indeed observed in 
Fig.lGJ where we show that Ke{ipii) has a "wavy" pattern 
at any fixed extraction radius. The scalar-wave nature 
of this pattern is well illustrated by the right panel of 
Fig. [6J There we take the time derivative of Re((/3ii) at 
the largest and smallest extraction radii, scaling the am- 
plitude of the signal by the ratio of the extraction radii 
(as expected for a wave scaling like 1/r), and we observe 
that the signal is indeed dipole scalar radiation emitted at 
twice the orbital frequency of the binary, consistent with 



the predictions of Sec. IIIB2 Furthermore, in Eq. (B15) 



of Appendix |B] we show that the amplitude of the ipu 
mode is consistent in order of magnitude with analytical 
predictions. We also observe a monopolc component (poo 
whose amplitude is comparable to the amplitude of ipn, 
consistently with analytical predictions. 

To understand the merger signal, when the two BHs 
collide and relax to a final nearly stationary state, it is 
useful to remember that, in vacuum, the merger of a 
BH binary with mass ratio q — 3 produces a Kerr BH 
with spin a/M ~ 0.543 |110j . Thus the lowest ringdown 
frequencies are expected to be, from perturbative cal- 
culations, Muj = 0.351 - 0.0936?, 0.476 - 0.0849i for a 
I = m = 1 scalar field, I = m = 2 gravitational mode 
[1111 I112| . We find that ^4 indeed rings down with 
Muj ^ 0.48 — O.OSli, in good agreement with pertur- 
bative calculations. An analysis of ipu yields a ringdown 
frequency Mtu ^ 0.36— 0.070i (with errors < 10%), which 
is roughly consistent with perturbative calculations of 
scalar (s = 0) perturbations of Kerr BHs |112| . Our sim- 
ulations also show that the ringdown phase, and indeed 
the entire scalar signal, scales with a. We conclude that 
our results indeed represent linear effects, as opposed to 
nonlinear mode couplings. 

Although not completely obvious, there is a small DC 
component in Fig. m (right panel) , which we estimate to 
be 






0.2, 



(77) 



at early times, where Iv'ii'* | is the absolute value of the 
waveform at a local peak (maximum or minimum) . This 
numerical estimate can be compared to the analytical 
prediction, Eqs. (B11)-(B12). We start by estimating 



the angular frequency x through the waveform frequency, 
and we find Mx ^ 0.025. Using Kepler's law, one can 
estimate the orbital separation, and these two ingredi- 
ents together with relations ( 32 ) allow us to estimate 



the relevant ratio of the time derivatives of expressions 
(B11)-( [bI2 ). We find a ratio which is smaller by almost 
one order of magnitude. This discrepancy can probably 
be explained by numerical uncertainties and strong-field 
nonlinear effects. 



It is apparent from Fig. [6] (left panel) that the Re{ipii) 
modes display a "drift" : after all the dynamics has died 
away, the field does not return to zero. Our data implies 

that at late times Re{ipii) 2.7x lO^'^r^i-^'' for Ma = 

2 X 10^^. The most natural interpretation of this drift is 
related to the DC component, Eq. (B5), which predicts a 
linear growth in time ~ roughly the same dependence that 
can be seen in Fig. [61 One should also bear in mind that 
the analytical result is a slow-motion expansion, whereas 
the numerical results cover only the highly dynamical, 
nonlinear merger signal; some deviation from a perfectly 
linear dependence is therefore expected. 

There are other possible contributions to such a drift. 
A second possible contribution is due to the nonvanishing 
of Im(i^ii) for the analytical solution (17) with 7 = 7r/2 
(recall that in both the single-BH solution (17) and in 
the numerical solution discussed here, the polar axis is 
orthogonal to the gradient, but the ordering of the axes 
is different, which explains why real and imaginary parts 
of the modes are swapped). However, the 1/r piece of 
that solution is orders of magnitude smaller than the am- 
plitude of the drift we observe numerically, and therefore 
unlikely to explain our observations. 

Another possible contribution comes from gravita- 
tional recoil. We are simulating an unequal mass binary, 
which acquires a kick from the origin of our coordinate 
axis as a result of the merger. The kick introduces "spu- 
rious" multipolar components with respect to a frame 
which is not comoving with the final BH. However, an 
ordcr-of-magnitude estimate shows that this effect is un- 
likely to explain the observed drift. In order of mag- 
nitude, the kick contribution to ip can be estimated by 

looking at the terms in Eq. (BIO), say 25rocoib that are 
proportional to vqm'- this yields 



V„ 



MaVcMVorh 

P/M 



(78) 



Here Ma = 2 x 10~^, the extraction radius is f/M ~ 10^, 
and forb is the orbital velocity. The maximum re- 
coil velocity from a nonspinning BH binary is of order 



WCM 



7 ■ lO^-* 1113J, so 



2?., 



< 10 ^^ even if v„ 



'recoil/ '^ 

approaches unity. This contribution to the dipole radia- 
tion is way too small to account for a significant portion 
of the drift seen in our simulations. 

Finally, a frame-dragging effect can also contribute 
with a nonzero drift for the scalar field. The coalescing 
binary drags the inertial frames, inducing a local rotation 
of the coordinate lines. This induces, near the binary, an 
apparent rotation in the y-z plane of the extracted scalar 
field, which determines a nonvanishing real part of (jjn. 
While the order of magnitude of this effect is roughly 
consistent with our numerical findings, the decay of the 
frame-dragging effect with extraction radius is not con- 
sistent with our data. Therefore frame dragging is not a 
dominant contribution to the observed drift. 

While the DC component accounts for the order of 
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magnitude of the drift observed in our numerical simu- 
lations, most likely the observed drift is due to a combi- 
nation of the effects mentioned above, and possibly oth- 
ers. In particular, by imposing constant-gradient bound- 
ary conditions at finite distance from the binary during 
the evolution we are effectively injecting energy into the 
system. This causes a growth of the scalar field, which 
may contribute significantly to the drift. This expecta- 
tion should be confirmed by longer simulations and/or 
by simulations where the "plates" generating the scalar 
gradient are located further away from the binary. We 
hope to return to this problem in future work. 



VI. CONCLUSIONS AND OUTLOOK 

We have investigated BH dynamics in external field 
profiles, by considering the very simple example of a 
constant scalar-field gradient. The broad features of our 
analysis should translate to other, more general settings, 
at least as long as the external force varies on length- or 
time-scales which are larger than the typical binary sepa- 
ration. Our results are in agreement with linear or slow- 
motion expansions, and show conclusively that black hole 
binaries evolving in a nontrivial background produce in- 
teresting scalar and gravitational-field dynamics. 

As discussed in Appendix |D] the scalar-field profiles 
currently considered in scalar-field dark matter models 
correspond roughly to Ma ^ 10^^" or less for a typical 
stellar-mass BH with M — 10 Mq. Because scalar ra- 
diation is proportional to the gradient, the experimental 
relevance of our setting for gravitational radiation from 
BH binaries, as observable by Advanced LIGO or similar 
instruments, seems negligible. However, it is interest- 
ing that gravitational-wave observations may yield up- 
per bounds on scalar field gradients at all. Furthermore, 
strong field gradients can be encountered in other - al- 
beit more speculative - dark matter configurations, such 
as supermassive boson stars [55] • Together with the fact 
that nature has an incredible knack to surprise us, the 
possibility to come across this type of signals should be 
seriously taken into account. 

Our analysis answers some questions and sparks many 
new ones: why exactly do large gradients develop insta- 
bilities? Are these instabilities of a physical or purely 
numerical nature? Even isolated BHs moving in a scalar 
field should accrete: how can we understand the details 
of this accretion process? Another interesting question 
concerns spinning BHs. The original linearized analysis 
by Press j95j shows that Kerr BHs should in principle 
align their rotation axis with the field gradient over long 
enough timescales. Numerical simulations of this align- 
ment and of its nonlinear development are an interesting 
(but numerically challenging) open problem, which prob- 
ably requires much longer simulations than those pre- 
sented in this work. 
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Appendix A: Boosted black-hole background 

In this Appendix we shall consider the solution of the 
Klein-Gordon equation in a boosted BH background, i.e., 
in the presence of a Schwarzschild BH moving with con- 
stant velocity v in the direction orthogonal to the charged 
planes. We shall show that this solution does not emit 
scalar radiation: only accelerated BHs moving in a uni- 
form scalar-field gradient can emit scalar radiation. 

To this aim, wc shall first consider the (regular) so- 
lution describing a scalar field on a Schwarzschild back- 
ground, generated by an infinite plane moving with veloc- 
ity V along the direction z orthogonal to the plane. Then, 
we shall boost back this solution, to obtain a moving BH 
and a scalar field generated by a fixed plane. 

The expression 



tp — 27r(T7 



{r - M) cos e + 2Mv 



V -r 
2M 



— loe 



2M 



(Al) 
describes a solution, regular at the horizon, of the Klein- 
Gordon equation (14) on a Schwarzschild background. 



Here r* = r+2AI \og(r/2M—l) is the tortoise coordinate, 
V ^ t + r^ is the standard advanced time coordinate, v 
is a velocity parameter and 7 = (1 — w'^)^^". For v ~ 
0, the previous expression reduces to the static solution 
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(16 1. For finite v and at asymptotically large distances Therefore 



it reduces to 



Lp r^ 2tt<t^{z + vt) , 



(A2) 



which is related to the asymptotic solution ( 16 1 by a sim 



pie boost. Therefore, Eq. (All describes a scalar field 



(on a Schwarzschild background) generated by an infi- 
nite plane moving with velocity —v along z. We shall 
now boost Eq. ( Al I and the Schwarzschild metric, in or- 



der to find a solution of the field equations describing 
a BH moving with velocity u in a field generated by an 
infinite charged plane at rest. The boost is more easily 
performed in isotropic coordinates (t,r,9,(f>), in which 



ip = 27rcr7 



-2Mv\og- 



1 



M/f + 



+ tv 
ir/M 



-4 + M/f + Af/M_ 
If we apply a Lorentz boost along the z-direction 

i = j{t + vz) 



(A3) 



we get 





z = -f[z + vt) 

X — X 

y = y, 


if = 27ro- 


z+\^{z-vt^ 

4 4- Mir + AtIM 




*= -4 + Mjr + 4f/M 



(A4) 



(A5) 



Note that the isotropic radial coordinate reads f^ = 
y^ +^'^{z-vlf-. 

Along the z— axis, this expression has the form 



'f 



2-Ka 



z + 



Af2(l-l6w) AP{l~lQv)vt 



Az Az"^ 

Introducing polar coordinates in the boosted frame 



0{-z-^). 
(A6) 





X = r sm a cos 




y — f sin 9 sin 




z — f cos 9 , 


we have (since 7^ 


- 1 = «2^2) 



(A7) 



f = \/r2(l + ^272 cos 61) - 2w72ircos6' + 72i;2t2 . (A8) 

If we assume small boosts v ^ 1, so we can neglect terms 
0{v'^), 7 ~ 1. We can also assume that vi <c f, i.e., 
that vi/f <^ 1, even tho ugh this quantity can be larger 



than V . Expanding Eq. ( A8 1 in these two dimensionless 



quantities, up to first order in w^ and second order in 
vt/f <^ 1), we find 



1 



vt 



cos 9 



(A9) 



if — 2TTcr 
AM^v 



- M^ 

r cos 9 -I I cos ( 

Ar 



2(^]cos^' 
r 



1-t- ( - \cos9 
r 



O 



1 



: 27rcr 

M'^vt 
' 2f2 



- M2 

r cos 9 H — (cos 9 — 16w) 

Ar 



cos9{cos9 — 8w) 



(AlO) 



As discussed in jST] j we find that the scalar charge is the 
(spherically symmetric component of the) coefficient of 
1/f in this expansion, divided by 2il/, thus 



Q 



-AnaMv . 



Furthermore, the scalar-field multipoles are 

M = SnaM^v , 

^ naAP 
V ^ — z. 



(All) 

(A12) 
(A13) 



All of these quantities are constant in time, therefore 
there is no emitted scalar radiation. We can conclude 
that if a BH moves with constant velocity in a scalar- 
field gradient, there is no scalar emission; in order to 
get scalar radiation the BH should have a nonvanishing 
acceleration, as shown analytically in Section |IIIB| and 
numerically in Section [VB[ 



Appendix B: Approximate solution for monopole 

and dipole radiation from a quasi-circular binary in 

a scalar gradient 

This Appendix completes the approximate solution 
discussed in Section |IIIB 2[ which describes the scalar 
field generated by a binary BH in quasicircular orbit in a 
scalar-field gradient. The center-of-mass acceleration is: 



OCM 



WCM 



^CM 



4a(°)x^°U g- 
(9 



Stq 



1) 



z cos X — y sin X 



(Bl) 



4a(0) (q - 1) 
l>^{q+l) 

4a(o) 



y{cosx — 1) + i 



(9-1) 



5rqx(*') {q + 1) 
+Vo{t~to), 



ysmx 



(0) 



smx 



(0) 



(B2) 



z(cosx^ ^ - 1) 



(B3) 



where q — M1/M2 is the mass ratio, x i^) = X ^+"00 is 
the zeroth-order orbital phase, which vanishes at t — t^^ 
and Vo — vcuit ~ to) is the initial velocity of the center 
of mass relative to the scalar gradient. Without loss of 
generality, the choice zcuit — to) — has been made. 
Note that in the equal-mass limit (9 — )• 1), the center-of- 
mass recoil vanishes, on account of symmetry. 
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Given the scalar charges 

_ STTCrAff 



Q2 = 



M 
M 



MvcM • 5 + M2 [d sin X + ax cos x] , 



MvcM ■ z — Ml [a sin x + ^X cos x] I , 



the monopole component of the scalar field is: 



M = Qi+Q2 = j^{ M{Ml + M|)(t/cM • z) + M^M^iMi - M2)[dsinx + ax cos x! 



while the dipole component reads 

V = QiZi{t) +Q2Z2{t) = X'cM +2?roLDC + 2?rcL c 

with 

SiraMilVhiMi - M2)a 



2?CM — 
2?rel, DC 



P, 



rcl, osc 



M 
SiraM^M^a 

M2 



zcMX cos X + (t^CM • ^) (y cos X + z sin x) 
az + axyj , 
d[ysin(2x) - icos(2x)] + ax[j/cos(2x) + isin(2x)] 



From the previous expression we find 



I? = I?CM+A-cl,DC+A- 



Cl, osc ; 



with 



2?CM = ^7 x(t'CM COS X - ^CMX Sin X) + (acM ■ z) [y cos x + z sm x) 



M 



^ 87rcrM?M|a, ._ ... „, , 
Aoi, DC = ^T9 I «^ + (2aX + aXjy , 



Aoi, 



Af2 

87rcrMfAf|a 
AP 



where M — Mi + M2 is the total mass, ucm and acM 
are the velocity and acceleration of the center of mass, 
respectively, and terms of order higher than l/c' in the 
dipole moment have been dropped. If we neglect radia- 
tive effects, and express the dipole waveform using the 
£, m multipole components introduced in (72), we find 



<^io = 0, (B13) 

ipi ±1 = zAe'F^^x _|_ gradient term , (B14) 

where the amplitude of oscillation is given by 



(B4) 

(B5) 

(B6) 

(B7) 
(B8) 

(B9) 
(BIO) 

(Bll) 
(B12) 



I 

and V = M1M2 1 M'^ is the symmetric mass ratio. For our 
simulation with q — ?, and Ma = 2 x 10~^ we find that 
Afx ^ 2 X 10^^, and therefore the theoretical prediction 
for the dipole amplitude is ^ ~ 4 x 10""'^°. This is in 
order-of-magnitude agreement with the observed dipole 
radiation in our numerical simulations. 



Appendix C: 3+1 Evolution equations 



+ (wcM ■ z)x{z cos X - y sin x) 



(d-2ax2)[ysin(2x) - zcos(2x)] + (4dx + ax)[y cos(2x) + zsin(2x)] , 



In terms of the BSSN variables defined in Eq. (651, 



5127r3. ^ M 2,,,.^2/3 

6 r 



(B15) 



the scalar field ip and the scalar curvature K^p defined in 
Eq. (44 1, the BSSN evolution equations are given by 
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dtX = rdmX + \x{o^K - d^r) , 



dtAj = rOmA, + 2i„(,a,)/3™ - ^Ajd,nr + X {aR^J - D^D.af' + a (.v ^,, - .^, ^„„ 



J^ An 2Ai An 



-ax { d^ip djip - -%-7™"9„i(^ d^^ \ 






X 



af^draK - ^"" ( ia"-^^ + 2a,„a ) - 3 (f * - 7™"f:„„) 5^/3'= - 4if^7™a™^. 



(CI) 
(C2) 

(C3) 

(C4) 

(C5) 



Likewise, we write the Hamiltonian and momentum 
constraints in terms of the BSSN variables as 



(3)i?^ + -K^ - ^^^^'A^^kAni = AKl + d.ipd'^ ,(C6) 



DmA"\ ~ '^d^K - |-i™.a,„x = 2K^d,^ . (C7) 



The lapse function and shift vector are evolved as in 
the case of vacuum general relativistic BH simulations 
(cf. |114| ) according to 






(C8) 
(C9) 



Following |115j . we use a position-dependent parameter 
77; specifically, we set 



i?2 



V^Vo^- 



n 



r2 



i?2 2(Mi|fi|+M2|f2|)' 



(CIO) 



where r is the coordinate distance from the origin, ri 2 
are the position vectors from either hole and Mi. 2 are the 
BH masses. Lapse and shift are initialized as a = ^/x ^^'^ 
/?' = 0, respectively. 



Appendix D: Order of magnitude of the scalar 
gradient in a cosmological scenario 

Scalar fields on galactic scales have been considered 
by many authors as a possible explanation for the rota- 
tion curves in galaxies and as alternatives to cold dark 
matter [52H5(Ij (see also [91j ) . The aim of this Appendix 
is to estimate the typical magnitudes of the scalar-field 
gradients predicted by these models. 

The rotational velocity near galactic centers Vg ^ 
200 km/s is approximately constant [116) . The mass 
M{r) contained within an orbit of radius r is roughly 



r 



given by 



M{t 



Then, in order of magnitude, the density is 



(Dl) 



(D2) 



on scales of the order r ~ 10 kpc = 3 x 10^^ km (see e.g. 
|90)). This density is provided by dark matter, or, in the 
scalar-field scenario, by the scalar field. 

One should be very cautious in comparing our station- 
ary, free scalar field configuration with those of cosmo- 
logical models. Indeed, in many of the works cited above 
the scalar field is rapidly oscillating, and the mass term 
and the potential always play a role. In any case, if our 
aim is to estimate the order of magnitude of the scalar- 
field gradient assuming that it yields galactic rotation 
curves, we can simply note that, neglecting the contribu- 
tion of the potential V{(j)) and restoring physical units, 
the mass-energy density is of the order (see e.g. |5F| ) 



Gp^\^,t\'+C^\^,r\\p^\<f>,t\'+C^\ 

therefore the gradient a ~ 4>^r is at most 



10"^^ km" 



10- 



10A/,r 



(D3) 



(D4) 



In our simulations the mass is normalized to the BH 
mass, i.e. we set M = 1. For an astrophysical BH {M = 
IQMq) moving near the galactic center a typical scalar- 
field gradient is therefore of the order Ma ^ 10~^°. This 
should be considered as a rough upper limit: in scenarios 
in which the scalar field is rapidly oscillating the kinetic 
term should contribute to the energy density more than 
the gradient. 
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